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1. INTRODUCTION 


The following is a brief survey of some recent developments in finite-element 
analysis technology which bear upon the three main research areas under consid- 
eration in this workshop: (1) analysis methods; (2) software testing and quality 
assurance; and (3) parallel processing. 

The variational principle incorporated in a finite-element computer program, 
together with a particular set of input data, determines the exact solution corre- 
sponding to that input data. Most finite-element analysis computer programs are 
based on the principle of virtual work. In the following we consider only programs 
based on the principle of virtual work and denote the exact displacement vector 
field corresponding to some specific set of input data by u EX • The exact solution 
u E x is independent of the design of the mesh or the choice of elements. Except for 
very simple problems, or specially constructed test problems, u EX is not known. 

We perform a finite-element analysis (or any other numerical analysis) because 
we wish to make conclusions concerning the response of a physical system to 
certain imposed conditions, as if u EX were known. We know the finite-element 
solution only which we denote by u FE . The solution u FE depends not only on the 
variational principle and the input data but also on the finite-element mesh and 
the choice of elements. We will assume that the finite-elements are exactly and 
minimally conforming and therefore the elements are completely characterized by 
their polynomial degree. We therefore control u FE by mesh design and the choice 
of the polynomial degree of elements. 

We wish to compute u FE so that u FE is close to u EX in some sense. For example, 
if we are interested in determining a stress intensity factor then we wish to have 
the stress intensity factor computed from u FE to be close to the stress intensity 
factor computed from u EX within some prespecified level of tolerance r. In general, 
we wish to determine functionals Vi(u FE ) (*' = l, 2...,n) so that: 


^»(vsx) ~ 
^■(vsx) 


< Ti 


(* = 1, 2, ... , n) 


( 1 ) 


The question naturally arises: how can we tell whether is close to $\(ubx) 

if we do not know u EX ? The answer is: by performing extensions. Both the 
estimation and control of error are based on extensions. 
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Extensions are systematic increases in the number of degrees of freedom either 
by mesh refinement, increase in the polynomial degree of elements or a combination 
of both. If the extension is by mesh refinement then the process is called h- 
extension* . If the extension is by increase in the polynomial degree of elements 
then the process is called p-extension**. If the extension is by a combination 
of proper mesh refinement and concurrent increase in the polynomial degree of 
elements then it is called h-p extension. Having performed an extension, we may 
draw conclusions concerning the overall quality of the approximate solution and 
the quality of any functional computed from u FB . 

2. OVERALL QUALITY 

The overall quality of approximation can be judged in terms of the estimated 
error in energy norm and errors in equilibrium. Estimation of error in energy 
norm is outlined in some detail and an example is presented. Procedures for 
assessment of the quality of approximation in terms of errors in equilibrium are 
briefly discussed. 

2.1. Estimation of error in energy norm. 

We know that the strain energy of the error U(u BX -u FB ) must decrease mono- 
tonically as we systematically refine the mesh or increase the polynomial degree 
of elements and a well developed, elaborate theoretical basis exists for the estima- 
tion of error in energy norm for the h-, p- and h-p extension processes. (See, for 
example, [1,2, 3, 4, 5].) The error in energy norm is defined as: 


II^EX - VFbIIb(O) = yU(^EX - «fs) (2) 

where fi represents the solution domain and U represents the strain energy. ||u £ x - 
«fe|U(o) is closely related to the root- mean-square of error in stresses [6], 

In the case of h- and p-extensions the estimate is of the form: 


||«£X - «Fe||e(0) < Jfp 


(3) 


where k and p are positive constants, N is the number of degrees of freedom. In 
the case of h-p extensions the estimate is of the form: 


\Gex-Zfe\\e ( o) < — 


(4) 


where k , and 9 are positive constants. These estimators are ’sharp’ for large N 
values hence the ’less than or equal’ (<) can be replaced by ’approximately equal’ 
(«) in (3), (4) when N is large. Therefore from (3) for large N values we have: 

log||usx - ufe||e(o) « log fc — /3 log N (5) 

If we plot log || u B x - «Fs|U(n) versus log N we see a downward sloping straight line. 
The absolute value of the slope is p, called the asymptotic rate of convergence. 
When p is large then the error decreases rapidly as A is increased. When p is 
small then the error decreases slowly. Of course, the error also depends on k which 


* h represents the size of elements. h-Extension involves letting h max —* 0. 

* p represents the polynomial degree of elements. p-Extension involves letting p m in —* 00 . 
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is generally not known a priori, but can be estimated from data obtained from 
properly performed extensions. This will be discussed later. When the estimate 
is of the form (3) the rate of convergence is said to be algebraic. 

When the estimate is of the form (4) and we plot log||u E x - UFs||s(n) versus 
log JV then for large N values we see a downward curving line [2, 3, 4, 5, 7]. In this 
case the rate of convergence is exponential: 

log || uex — ufe ||s(0) « log k — 7 (log e) N $ (6) 

where e is the base of the natural logarithm. If we plot log||u B x - «FBlU(n) versus 
N e (not log N as before) then we see a downward sloping straight line. It is known 
that under conditions which are generally satisfied in practice 6 > 1/3 [5]. 

All error estimation techniques are based on extension. Because in general 
the exact solution u EX is not known, the only information available to us is how 
the finite-element solution u FE behaves when the number of degrees of freedom is 
increased either through mesh refinement or increase in the polynomial degree of 
elements. Such information, together with an estimate or hypothesis concerning 
the magnitude of the error, or its rate of change with respect to N, is essential to 
all error estimation. Of course, the estimate or hypothesis must be asymptotically 
correct: as N — ► oo the estimated error must approach zero at the same rate as the 
true error does. Therefore the quality of error estimators should increase with N. 

P-extension makes it convenient and inexpensive to obtain information con- 
cerning the rate of change of U(u FE ) with respect to N. In the p-version hierarchic 
basis functions are used. Therefore the stiffness matrices and load vectors corre- 
sponding to polynomial degree p are embedded in the stiffness matrices and load 
vectors of polynomial degree p + 1. Once a solution is available for polynomial 
degree p max , all solutions corresponding to p = 1, 2,...,p max - 1 can be readily and 
inexpensively obtained. Specifically, we write: 


|| UEX - U FE \\%(ft) = U(u E x - UFe) = | U(u E x) - ^(«Fs)| « 
Let us assume for the moment that U(u EX ) > U(u FE ). In that case: 

U{u EX ) - U(u FE ) » 


(7) 

( 8 ) 


We have three unknowns: U(u E x), k and /?. If we have three values of U(u FE ) and 
N corresponding to three different values of p, then we have three equations for 
computing the unknowns. Let us denote these three values by U p , U p - 1, U p - 2 and 
N p , JVp_ i , Ap_ 2 and U(u EX ) by U. Then from (8) we have: 


log 


log 


V^Up 
U-U p _! 
U-U P - 1 
u - E/p_ 2 




log 


log 


Np-i 

A P 
Np - 2 

Np—i 


(9) 


Denoting the right hand side of (9) by Q, we have: 


U-U p 
U - U p . i 




( U-Up- 1 \ 

\U-U P . 2 J 


Q 


( 10 ) 


To obtain an estimate of the exact strain energy U, we need to solve (10). The 
solution is expected in the neighborhood of U p . Because convergence of the strain 
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energy is monotonic, we know that U > U p when U p > U p - 1 . Conversely, U < U p 
when U p < U p _ 1 . Eq.(lO) would not be different if U(u EX ) < U{u FE )\ therefore, the 
restriction that U(u EX ) > U(u FE ) is not essential. Computational experience has 
shown this estimate to be reliable and generally accurate, with the accuracy of the 
estimate increasing with the accuracy of U p . 

2.2. Example. 

The following test problem is representative of plate and shell intersections 
and reentrant corner problems in general. An L-shaped plane elastic body of 
thickness t is loaded by tractions. The tractions are computed from a stress field 
which satisfies the equilibrium and compatibility equations and the stress free 
conditions along the reentrant edges. Specifically, the stress field corresponds to 
the first (symmetric or ’Mode 1’) term of the asymptotic expansion of u EX about 
the reentrant corner. (See, for example, [8].) Therefore the exact solution is 
known. Specifically, the components of u EX in the coordinate system shown in 
Fig. 1 are: 


u x = r x [(/c — Q(X +1)) cos A 0 — Acos(A — 2 ) 0] ( llcz) 

2 G 

u y ~ r * [( K + + ^)) sm A 0 + A sin(A — 2) 6] (H&) 

2 G 

where A is a generalized stress intensity factor; A = 0.544483737; Q = 0.543075579; G is 
the modulus of rigidity and k depends on Poisson’s ratio u only. For plane strain: 
K = 3 - 4 1 /. We assume plane strain conditions and v = 0.3, therefore in this case 

K = 1 . 8 . 



Fig. 1. L-shaped plane elastic body. 


The stress tensor components are: 

a x = A A r A_1 [(2 - Q(A + 1)) cos(A - l) 6 - (A - 1 ) cos(A - 3) 0] (12a) 

<jy = A A r A_1 [(2 + Q(A + 1 )) cos(A - 1 ) 6 + (A - 1 ) cos(A - 3) 0] (126) 

r xy = A A r A_1 [(A — 1 ) sin(A — 3) 0 + Q(X + l) sin(A — 1 ) 0] . ( 12 c) 

i 
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Because we know the exact displacement and stress fields we can compute the 
strain energy of the exact solution: 

A 2 a 2X t 

U(uex) — 4 . 15454423 — — — ( 13 ) 

E 

where E is the modulus of elasticity. The relative error in energy norm is defined 
as follows: 



Using the mesh shown in Fig. 2 finite-element solutions were obtained for 
p=l to 8. The computations were performed by a new computer program, called 
PROBE [9], The number of degrees of freedom, the computed strain energy, the 
estimated and true relative errors in energy norm, computed from eq. (14), are 
shown in Table 1. 


22 



Fig. 2. Mesh design. 


The results presented in Table 1 are typical of the quality of the error estimate 
we can obtain by means of the procedure described in Section 2.1. When the mesh 
is strongly graded toward the point of singularity then the convergence path (the 
log(e r ) e versus logiV curve) looks like an inverted S [3,4,5,10]. For low N values 
the rate of convergence is nearly exponential and the downward slope increases 
with N. In this segment the estimated error is conservative. Near the inflection 
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Table 1. Estimated and true relative error in energy norm. 


p 

N 

U(ufe)E 
A 2 a 2 * t 

2/3 

Est.'d 

( e r)B 

True 

(er)fi 

1 

41 

3.8860880 

— 

— 

25.42 

2 

119 

4.1248326 

— 

— 

8.46 

3 

209 

4.1481150 

1.93 

5.34 

3.93 

4 

335 

4.1526504 

2.76 

2.02 

2.14 

5 

497 

4.1536354 

3.05 

1.01 

1.48 

6 

695 

4.1539746 

2.45 

0.80 

1.17 

7 

929 

4.1541390 

1.83 

0.75 

0.99 

8 

1199 

4.1542378 

1.39 

0.75 

0.86 

oo 

oo 

4.1545442 

1.09 



0 


point, (i.e. where the curvature of the convergence path changes from negative to 
positive) the estimate is the least accurate and not conservative, nevertheless as 
we see in this example, it remains close. The estimate then becomes progressively 
more accurate as the asymptotic range of the p-extension is entered. In this case 
the correct asymptotic rate of convergence is /? = A = 0.5445. At p = 8 the computed 
value of fi is approximately 0.7 with (3 decreasing. 

2.3. Equilibrium tests. 

Smallness of error in energy norm is a necessary but not sufficient condition 
for ensuring that the overall quality of the finite-element solution is good. It is 
possible to produce examples where the estimated error in energy norm is small 
(under 1 percent) yet the error in overall equilibrium is large (well over 10 percent) . 

Although we do not know u EX we know that u EX satisfies the equations of 
equilibrium and the law of action and reaction. We can, therefore, assess the 
quality of the finite-element solution by examining to what degree u FE satisfies 
equilibrium and the law of action and reaction. Specifically, we can perform: (l) 
overall equilibrium tests; (2) element by element equilibrium tests and (3) action- 
reaction tests. 

In the overall equilibrium test we ’cut’ the structure from its supports and 
integrate the tractions, computed from the u FE , to obtain the reactions. In this 
way a free body diagram is produced. The error in equilibrium must be small in 
relation to the magnitude of applied forces. For example, we can define: 


F = 


N 



(15) 


where r, are the (global) load vector components. The term F is a suitable measure 
of the magnitude of the applied load. 

In the element-by-element equilibrium test individual elements (or any group 
of elements) are separated from the model and tested for equilibrium. Specifically, 
we denote the element domain by n e and its boundary by dfl e . We compute: 


qj e * = f (crijj + X { ) tdx 1 dx 2 + I tJijrijtda (i,j = 1, 2) 
J n, Jdn e 


(16) 
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where are the stress components computed from the finite element solution; X { 
represents the applied body force components and n } represents the unit normal 

to dd e . The summation convention is used. We should have |gj e) | = \J q\ e W^ as 
well as the absolute values of each component of both expressions on the right 
hand side of (16) small in relation to the magnitude of the applied loads. If |gj e) | is 
small but the absolute values of the integral expressions are not small then there 
is a local error but, according to Saint-Venant’s principle, the effect of the local 
error will not be substantial at some distance from the element in question. If 
|g( e )| is large, even after p-extension was performed, then element e, and possibly 
its neighbors, should be subdivided. Thus the element by element equilibrium 
test provides information about the quality of mesh design. In many cases minor 
local refinement (for example, dividing one element into two elements) can have a 
highly beneficial effect on the overall quality of approximation when p-extension 
is used. 

In the action-reaction test we compute the stress resultants along interelement 
boundaries and external element boundaries where tractions are applied. Along 
interelement boundaries the stress resultants computed for neighboring elements 
should have nearly the same absolute value and opposite sense. Along external 
boundaries the resultants of the applied tractions and the tractions computed from 
the finite element solution should be nearly the same. 

Examples of equilibrium tests are presented in [11]. 

3. LOCAL QUALITY 

Having ascertained that the overall solution quality is acceptable, we are ready 
to compute the quantities which are of principal interest, i.e. Smallness of 

error in energy and equilibrium does not guarantee that all functionals computed 
from u F e are accurate. It is advisable to perform convergence tests on at least 
the more difficult functionals. We demonstrate the procedure by computing the 
direction and magnitude of the principal stresses at a point close to the reentrant 
corner. We selected the point r = 0.025a; 6 = 30° in the coordinate system shown in 
Fig. 1. The stress components computed from the exact solution are: 

o x = 3.67198i4a A_1 a y = 7.68696 Aa A_1 r xy = 0.698375Aa A_1 (17) 

Therefore the principal stresses <xi and <r 2 and the direction of the first principal 
stress ax from the positive x-axis, denoted by are: 

cr ! = 7.804 Ao a_1 <7 2 = 3.554i4a A_1 6i = 80.4° (18) 


In general functionals, other than the strain energy, do not converge mono- 
tonically, nevertheless the fact that convergence has occurred should be obvious. 
Here cr x and cr 2 happen to converge monotonically but does not. We see that the 
state of stress is known with sufficient accuracy for engineering purposes at p=4 
(335 degrees of freedom, see Table 1). Extension beyond p=4 merely confirms 
that convergence has occurred to within the range of precision normally expected 
in engineering computations and thereby establishes reliability of the data. 

This test problem demonstrates that accurate stress data can be obtained in 
the very close proximity of stress singularities. Other examples and additional 
discussion of this point are presented in [10,12]. 
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Table 2. Principal stresses at r = 0.025a; 6 = 30°. 



cri 

CT2 

0i 

p 

Aa A-1 

Aa x_1 

(degrees) 

1 

7.096 

2.854 

92.1 

2 

7.441 

3.047 

84.1 

3 

7.532 

3.294 

81.6 

4 

7.731 

3.483 

80.9 

5 

7.754 

3.525 

80.6 

6 

7.773 

3.545 

80.5 

7 

7.786 

3.551 

80.5 

8 

7.791 

3.553 

80.4 

oo 

7.804 

3.554 

80.4 


It is possible also to compute various functionals from u FE using advanced 

methods of extraction [13,14,15] -For example, we may wish to determine the value 

of the generalized stress intensity factor A [see eq.’s (11a, b), (12a, b,c)]. Such 

procedures based on [13,14,15] have been implemented in PROBE [9]. 

4. CONCLUSIONS AND RECOMMENDATIONS 

(1) Extensions are essential for both the estimation and control of error in finite- 
element computations. 

(2) We are in a much better position today than we were, even just one year 
ago, from the point of view of understanding how an advanced finite element 
software system should be designed so that (a) the solution is obtained at very 
nearly the theoretically optimal efficiency and (b) the user is provided with the 
capability to estimate and control the quality of engineering data computed 
from the finite-element solution at a small marginal cost. This is because now 
we understand the interplay between mesh design and the polynomial degree 
of elements. 

(3) P-extension, coupled with properly graded meshes, is the most efficient method 
for controlling error in finite-element computations. 

(4) The proper mesh design is such that points of singularity (and areas where the 
solution changes rapidly over short distances) are isolated by one or more layers 
of small elements, with the elements graded in geometric progression toward 
the points of singularity. In this way both the global and local behavior of the 
solution can be represented without compromising the accuracy of either. 

(5) Implementation of advanced extraction methods for the the computation of 
certain engineering data, such as stress intensity factors, will further increase 
the efficiency and reliability of computations. 

(6) The p-version is well suited for implementation on parallel processors because 
the data are organized in relatively few, large units. This logical organization 
reduces the overhead associated with parallel processing. 

(7) The substantial increase of efficiency in finite-element computations through 
the use of h-p extension and the availability of parallel and vector processing 
technology make it possible and desirable to model plate and shell problems 
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using hierarchic sequences of plate and shell theories in conjunction with fully 
three-dimensional representation. The various plate and shell theories are 
nothing more than specializations of the three-dimensional theory of elastic- 
ity through restrictions imposed on the variation of the displacement field in 
the direction of the normal. Such restrictions generally do not hold near sup- 
ports, stiffeners, cut-outs, plate and shell intersections, etc. which are the 
areas where cracking and delaminations originiate and therefore of the great- 
est concern to analysts and designers. These areas can be properly modeled 
by three-dimensional representation only. The use of hierarchic extensions to- 
ward higher order plate and shell theories will permit us to assess and control 
the quality of approximation in relation to three-dimensional theory. 

(8) Although linear theory is properly the first and most generally used approach 
to structural modeling, it should be possible to ascertain by a posteriori anal- 
ysis whether engineering conclusions drawn from a numerical model would be 
different if geometric and material nonlinearities were considered. We can view 
linear theory as the simplest of a hierarchic system of theories. Much the same 
way as we estimate error by the use of extension processes within the frame- 
work of linear theory, we should be able to estimate error by extension within 
the hierarchic system of theories. This important area has not received much 
attention in the past. Because it bears on the reliability of computed data, and 
the engineering conclusions based on them, it deserves serious consideration. 

(9) In some areas our ability to compute data is already greater than the material 
scientists’ ability to tell us what data should be computed. For example, it 
is not fully understood what parameters govern crack initiation. The reason, 
at least in part, is that the conventional finite-element method tends to yield 
’fuzzy’ data in areas where stresses change substantially over short distances. 
Proper use of h-p extension, coupled with advanced extraction methods, per- 
mits us to compute any stress field parameter with arbitrary precision. This 
removes an uncertainty from the phenomenological characterization of mate- 
rial response to various stress fields. Of course, such characterization can be 
developed only through joint experimental-analytical investigations. 
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